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Abstract 

We describe a new algorithm for computing the Voronoi diagram of a set of n points in 
constant-dimensional Euclidean space. The running time of our algorithm is (9(/lognlog A) 
where / is the output complexity of the Voronoi diagram and A is the spread of the input, the 
ratio of largest to smallest pairwise distances. Despite the simplicity of the algorithm and its 
analysis, it improves on the state of the art for all inputs with polynomial spread and near-linear 
output size. The key idea is to first build the Voronoi diagram of a superset of the input points 
using ideas from Voronoi refinement mesh generation. Then, the extra points are removed in a 
straightforward way that allows the total work to be bounded in terms of the output complexity, 
yielding the output sensitive bound. The removal only involves local flips and is inspired by 
kinetic data structures. 



1 Introduction 

Voronoi diagrams and their duals, Delaunay triangulations, are ubiquitous in computational geom- 
etry, both as a source of interesting theory and a tool for applications [5J. Algorithms for planar 
Voronoi diagrams abound where optimal algorithms are known. In higher dimensions, the situa- 
tion is complicated by the large gap between the best-case and worst-case output complexity. A 
Voronoi diagram of n points in M. d can have between 0(n) and Q(n^ d ^) faces |2Uj (we suppress 
constant factors that only depend on d). This motivates the search for output-sensitive algorithms 
and analysis, where the time and space guarantees depend on both the input size n and the number 
of output faces /. 

Computing Voronoi diagrams in M rf reduces to computing convex hulls in This relation- 

ship is perhaps most clear in the dual, where the Delaunay triangulation is the projection of the 
lower hull of the input points lifted into M. d+1 by the standard parabolic lifting: 



All of the previous literature on output-sensitive Voronoi diagrams in higher dimensions is directed 
at solving the more general problem of computing convex hulls 

Seidel gave an algorithm based on polytope shelling that runs in 0(n 2 + f log n) time [Jjj]. The 
quadratic term is from a preprocess that solves a d-dimensional linear program for each input point. 
Matousek and Schwartzkopf showed how to exploit the common structure in these linear programs 
to improve the running time to (^(n 2-2 /^/ 2 ^ 1 ) log ^ n + / logn) 

Another paradigm of algorithms uses the dual notions of gift- wrapping [22J and pivoting [4j. 
Both approaches can enumerate the facets of a simple polytope in 0(nf) time, thus paying ap- 
proximately linear time per face. Since Voronoi diagrams have at least n faces, these methods do 
not improve on the LP-based methods except that they avoid the exponential dependence on the 
dimension inherent in such methods. 

Chan gave an algorithm based on gift- wrapping that runs in 0{n log /+(n/) 1— 1 /^ d / 2 l+ 1 ) log ^ 1 ) n) 
time [6]. A later work by Chan et al. gave an algorithm that runs in 0((n + (n/) 1 - 1 /r(^+i)/2l + 
j n i-2/\(d+i)/2\^ logCK 1 ) n ) time [7J. Even when / = G(n), the running time is poly(n) per face. 

Better bounds are known for 3- and 4-dimensional Voronoi diagrams where truly polylogarithmic 
time per face algorithms are known. Chan et al. gave an algorithm that achieves 0(/log 2 n) for 
M 3 [7]. Amato and Ramos gave an 0(f log 3 n)-time algorithm in IR 4 [2\. 

Voronoi diagrams and Delaunay triangulations are used in mesh generation (see the recent 
book by Chan et al. [10J). Extra vertices called Steiner points are added in a way that keeps the 
complexity down. Perhaps surprisingly, the number of vertices increases, but the total number 
of faces can decrease. Such a mesh can be constructed in 0{n log A) time using only 0(n log A) 
vertices [131 US] , 

In this paper, we propose a new algorithm for constructing Voronoi diagrams that uses Voronoi 
refinement mesh generation as a preprocess. Then, it removes all of the Steiner points using 
a method derived from the field of kinetic data structures. We prove that at most 0(/logA) 
local changes occur during the removal process. Each local change requires only constant time to 
process. These local changes are ordered via a heap data structure which adds an extra factor 
of 0(log(/log A)) = 0(logn + log log A) to the running time. We assume an asymptotic floating 
point model of computation where points are represented by their coordinates with 0(logn)-bit 

To avoid confusion, when reporting running times of known algorithms, we always give the results as they apply 
to Voronoi diagrams rather than convex hulls. 
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floating point numbers [13J . In such a model, log log A = O(logn). Thus, the total running time is 
0(/ lognlog A). 

Unlike previous work on output-sensitive Voronoi diagram construction, our algorithm does 
not use a reduction to the convex hull problem. Instead, it uses specific properties of the Voronoi 
diagram to get an improvement. 

Contribution We present the MeshVoronoi algorithm, a new, output-sensitive algorithm for 
computing Voronoi diagrams in M d . MeshVoronoi runs in 0(/ log n log A) time. When A = 
poly(n) such as the case when input points have integer coordinates with O(logn) bits of precision, 
the running time is 0(/log 2 n). Even the 0(n^ d ^ 2 ^) worst-case inputs for Voronoi diagrams can be 
represented with polynomial (even linear) spread. We get an improvement over existing algorithms 
for inputs with polynomial spread in dimension d > 3 when / = o(n 2 ~ 2 /l" d / 2 l ). Moreover, the 
analysis in terms of the spread is often quite loose. The running time really depends on the aspect 
ratios of the individual Voronoi cells of the output, for which the spread stands in as an easy to 
describe upper bound. The other advantage of the MeshVoronoi algorithm is that it is simple 
both to describe and to analyze. 

Related Work We will make use of the flip-based construction of weighted Delaunay triangu- 
lations similar to that presented by Edelsbrunner and Shah [TT]. In that paper, the concern was 
to add a single point to a regular triangulation but when run backwards, it describes the removal 
of a single point by local flips. We are interested in removing sets of points simultaneously, the 
Steiner points of the mesh. In this respect, the problem more resembles the Delaunay triangulation 
splitting problem for which Chazelle et al. gave a linear time algorithm for the plane The only 
higher dimensional analogue of this result is the extension by Chazelle and Mulzer to the case of 
splitting convex polytopes in M 3 |9j. 

Our algorithm may be viewed as a special case of a kinetic convex hull problem, and indeed, 
the main tools come directly from the literature on kinetic data structures [12J. In general, the 
kinetic convex hull problem is much harder than the instances arising in our algorithm and it is 
only because of the specific geometric structure of these instances that we are able to prove useful 
bounds. 

Joswig and Ziegler presented a different approach to output sensitive convex hulls using homol- 
ogy calculations [T7]. This algorithm is shown to be output-sensitive for simplicial polytopes like 
those produced for Delaunay triangulations of points in general position. However, they do not 
improve bounds on the asymptotic running times in terms of n and / as their construction passes 
through several reductions. 

Many low-dimensional algorithms for computing Voronoi diagrams depend on incremental con- 
struction, where the points are added one at a time. When the input can be degenerate, Bremner 
showed that incremental constructions cannot be output-sensitive [5] . Our algorithm does resemble 
an incremental algorithm. The main difference is that it adds more than just the input points in 
the first phase of the algorithm. The harder work is removing these extra points 

2 Background 

Points and Distances in Euclidean Space We will deal exclusively with the case of points in 
d-dimensional Euclidean space. The Euclidean norm of x is denoted ||x|| and thus the Euclidean 
distance between two points x and y is \\x — y\\. For a point x E M and a compact set U C M. d , 
define d(x, U) := mm ye u \\x — y\\. 
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The spread A of a set of points is the ratio of the largest to smallest pairwise distances. We 
will assume an asymptotic floating point notation so that coordinates are stored as floating point 
numbers with 0(logn)-bit words (see [13] for a more detailed treatment of this model). In such a 
model, the spread is 2°^ in the worst case. 

Voronoi Diagrams and Power Diagrams Let P C l d be a finite set of points. The Voronoi 
diagram of P is a cell complex decomposing M. d into convex, polyhedral cells such that all points in 
a cell share a common set of nearest neighbors among the point of P. Formally, for a subset S C P, 
the Voronoi cell Vorp(5) is defined as the set of all points x of M. d such that d(x, P) = \\x — y\\ for 
all y G S, i.e. 

Vorp(S) := {x G R d : d{x, P) < \\x - y\\ for all y G 5} 

When S = {v} is a singleton, we will abuse notation and refer to its Voronoi cell as Vorp(?;) instead 
of Vor For points in general position (no k + 3 points lying on a common fc-sphere), the 
affine dimension of Vorp(S') = \S — 1|. 




Figure 1: A weighted Voronoi diagram with weights indicated by disks. As the weights of some 
points increase from left to right, some cells disappear from the diagram altogether. 

The Voronoi diagram of P, denoted Vorp, has a natural dual called the Delaunay diagram, 
denoted Delp. For point sets in general position, the Delaunay diagram is an embedded simplicial 
complex called the Delaunay triangulation. The Voronoi/Delaunay duality is realize combina- 
torially by inverting the posets of the corresponding cell complexes, identifying each /c-face of the 
Voronoi diagram with a (d — fc)-simplex of the Delaunay triangulation. 

If the points P have real-valued weights w : P — > R and the Euclidean distance from p G P to 
x G M. d is replaced with the power distance ir p {x) 2 = \\x — p\\ 2 — w(p) 2 , then the Voronoi diagram 
becomes a weighted Voronoi diagram, also known as a power diagram. The dual is still 
well-defined and is known as the weighted Delaunay diagram (or the weighted Delaunay 
triangulation when it is a simplicial complex). 

Weighted Delaunay diagrams are the projections of convex polytopes in one dimension higher. 
In fact, interpreting the power distance to the origin as a height function for the points lifted into 
gives the weighted Delaunay diagram as the projection of the lower convex hull of these lifted 
points. In particular, setting all weights to gives the (unweighted) Delaunay diagram as a convex 
hull in Thus, there is a strong connection between the problems of computing convex hulls, 

Delaunay triangulations, and Voronoi diagrams. 

Flips in Triangulations Bistellar flips are a useful way to make local changes in triangulations. 
This operation takes a collection of d + 2 vertices S whose convex hull ch(5) is a subcomplex of the 
triangulation, and replace its interior with a new triangulation. In M 2 , there are three classes of 
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flips: those that swap the diagonal of a quadrilateral, those that introduce a new vertex in a triangle 
by splitting it into three, and those that remove a vertex incident to exactly three triangles. These 
are called respectively (2, 2)-, (1, 3)-, and (3, l)-flips, where the numbers correspond to the number 
of d-simplices removed and inserted respectively. In d dimensions, there are similar (i,j)-fhps for 
nonnegative integers i, j such that i + j = d + 2, where each replaces i ci-simplices with j new 
d-simplices. 

Any pair of weighted Delaunay triangulations can be transformed, from one to the other by 
a sequence of bistellar flips. This process of local changes is at the heart of incremental con- 
structions Any continuous change in the weights of a set of points causes a change in the 
corresponding weighted Delaunay triangulation that can be realized by a sequence of bistellar flips. 
In such a case, the flips happen exactly at those moments when the diagram is no longer a tri- 
angulation. For unweighted points, this happens when some set of d + 2 points lie on a common 
(d — l)-sphere. This can be tested by a linear predicate. That is, d + 2 points pi, . . . ,Pd+2 he on a 
(i-sphere if and only if 
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These represent the degenerate configurations of points. As before, for weighted points, we replace 
the norm with the power distance to find that d+2 weighted points form a degenerate configuration 
if and only if 

Pi ■■■ Pd+2 

\\Pi\\ 2 ~ w{pif ■■■ \\p d+2 \\ 2 - w{p d+2 f 

1 ••• 1 
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Voronoi Aspect Ratios and Voronoi Refinement The in-radius of a Voronoi cell Vorp(g) 
is the radius of the largest ball centered at q contained in Voi p(q). The out-radius of Vorp(g) 
is the radius of the smallest ball centered at q that contains all of the vertices of Vorp(g). Such a 
ball contains all of Vorp(g) for bounded Voronoi cells. The aspect ratio of the Voronoi cell of q is 
the ratio the out-radius over the in-radius, denoted aspectp(g). A Voronoi diagram has bounded 
aspect ratio if every cell has bounded aspect ratio. More generally, we say a set of points M is 
r-well-spaced if for all v G M, aspect m(v) < r. 

Bounded aspect ratio Voronoi diagrams have many nice properties. The most relevant for our 
purposes is that no ci-dimensional Voronoi cell has more than 2°^ facets (faces of codimension 1). 

For any set of n points P, there exists a r-well-spaced superset M for any constant r > 2. 
Moreover, such a superset can be computed in 0(nlogn + |M|) time [18] with \M\ = 0(log A) [15] . 
The process of adding points to improve the Voronoi aspect ratio is called Voronoi refinement. 
It is perhaps more widely known in its dual form, Delaunay refinement. The extra points added 
are called Steiner points. 

The analysis of Voronoi refinement depends on a function called the Ruppert feature size, 
defined for all x G R d as the distance to the second nearest input point. 

fp(x) := maxd(x, P \ {p}) 
p&p 

One defines the feature size with respect to a set M similarly and denote it ¥m ■ Voronoi refinement 
produces meshes r-well-spaced points M such that for each vertex v G M, 

fpO) < Kf M (v), 
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where K = . 

The total number of output points is, up to constant factors, determined by the feature size 
measure of the domain 0, defined as 

m := / 

Jn ipW 

For a wide class of input domains, the feature size measure is 0(n) [21J. For general inputs, the 
bound of 0(n log A) mentioned above is well known and can even be derived as a corollary to 
Lemma [5] below. 



Sparse Voronoi Refinement The meshing preprocess that we use is called Sparse Voronoi 
Refinement (SVR) and is due to Hudson et al. [H]. SVR is able to avoid the worst case complexity 
of Voronoi diagrams because it guarantees that the intermediate state is always a well-spaced point 
set and thus has a Voronoi diagram of linear complexity. 

For the case of point sets in a bounding box, the SVR algorithm is easy to describe. It is an 
incremental construction that proceeds by alternating between two phases called break and clean. 
The break phase attempts to add a Steiner point q at the farthest vertex of a Voronoi cell that 
contains an input point that has not yet been inserted. If there is an input point p too near to q, 
then p is added instead. The clean phase repeatedly attempts to add the farthest vertex of any cell 
with aspect ratio greater than r until none are left. As in the break phase, if ever there is an input 
point too close, then it is added instead. 

The running time of SVR is 0{n log A + \ M\). The nlog A term comes from the point location 
data structure which associates each point with the Voronoi cell that contains it. When attempting 
to add a point, the nearby Voronoi cells are checked to see if any input points are nearby to insert 
instead. 

Acar et al. developed an efficient implementation of SVR in 3-d pQ. It can also be efficiently 
parallelized [15]. Recently, Miller et al. showed that a variation of SVR runs in 0(n log n + \M\) 
time by using a more complex point location scheme [18] . 



3 Algorithm 

In this section, we describe the MeshVoronoi algorithm. It has three phases: a preprocess where 
the input points are placed in a bounding box and meshed, a removal phase that eliminates most 
of the Steiner points by incremental flipping, and a cleanup phase that removes the outer bounding 
box. The main data structure is a heap that stores the facets that are scheduled for removal by 
flipping. We call this the flip heap. 

Throughout the algorithm, there is a global time parameter t. There will be a superset M of 
the input points P. At time i, Mf denotes the set M with weights, where the weight of a point p 
at time t is given as 

w( t)-={ ^ ifpG P 
' 1 otherwise 

Using this weighting scheme, there exists a sufficiently large t such that for all t' > t, Debvfi = 
DelA/ t , . We use M* to refer to Mt, where t is some such sufficiently large value. 



3.1 Checking Potential Flips 

A set S of d + 2 points forms a potential flip at time t if ch(S') C Deljv/i- A potential flip 
can be identified with any of its interior facets. For example, in the plane, (2, 2)-flips are usually 
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Figure 2: An illustration of the algorithm from left to right. Starting with a point set, it is 
extended to a well-spaced superset. The cells of the Steiner points are shaded. Then, the weights 
of the input points are increased until the extra cells disappear. 



associated with the edge in the interior of a convex quadrilateral that will get replaced during the 
flip. In general, the facet representing a potential flip stands for the d + 2 points comprising the 
two simplices sharing that facet. 

The weights of the points of M vary with with the time parameter t. For p £ M, the power 
distance from x 6 M. d to p at time t is given as 
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if p G P 
otherwise 



So, for d + 2 points pi, ■ 
value of t such that 
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The determinant on the left hand side is a linear function of t, so it is easy to solve for t. When 
evaluating a potential flip, the algorithm simply checks if the time associated with the flip is before 
or after the current time. This approach to viewing linear, geometric predicates as polynomial 
functions of time is well established in the area of kinetic data structures [12]. In our case, only 
one row of the matrix is changing with time, and the change is linear, so there is no need to solve 
higher degree polynomial systems as is common in other kinetic data structures problems. 



3.2 Preprocessing 

The first step in the algorithm is to add a constant number of points to form a bounding region 
around the input. This can be done so that the total complexity of the convex hull of the augmented 
point set is a constant. Second, the Sparse Voronoi Refinement algorithm adds 0(nlog A) Steiner 
points to produce a r- well-spaced superset M, for a constant r > 2 (choosing r = 3 is reasonable). 
Third, the facets of the Delaunay triangulation of M are each checked for a potential flip and added 
to the flip heap accordingly. 



3.3 Flipping out Steiner points 

We maintain the weighted Delaunay triangulation through a sequence of incremental flips induced 
by the changes in weights. While the flip heap is nonempty, we pop the next potential flip. The time 
t is set to be the time of this potential flip. If the flip is still valid, i.e. all of the relevant facets are 
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still present in the weighted Delaunay triangulation at time t, then we perform the flip. Otherwise, 
we do nothing and continue. If any new facets are introduced, we check them for potential flips 
and add them to the flip heap accordingly. Then we loop. 

3.4 Postprocessing 

To complete the construction, the algorithm removes all boundary vertices and all incident Delaunay 
simplices. 

4 Analysis 

4.1 Boundary Issues 

The weighting will not remove all of the Steiner points. However, as the following lemma shows, 
only the boundary vertices will remain. 

Lemma 1 (Only boundary vertices remain). If q G Del^t for some Steiner point q and all t > 0, 
then q is a boundary vertex. 

Proof. Let p be any input point and let q be any non-boundary Steiner point. Let 

t* = max \\p — x\\ 2 . 

rreVor M (<?) 

Such a maximum exists because the Voronoi cells of non-boundary vertices are compact. Suppose 
for contradiction that q G DelM t for some t > i*. Since q G DeiM t , there must exist some point 
y G YoTM t (q)- It follows that y G VorM(<z) as well and so \\p — y\\ 2 < t±. We now observe that 

Mv> tf = \\P ~ 2/11 2 - *<**-*<0<||9- y\\ 2 = n q (y, t) 2 , 

contradicting the assumption that y G VorM t (q) • □ 

We need to show that removing the boundary vertices in a naive way as a post-process gives 
the correct output. For this, it will suffice to show the following lemma. 

Lemma 2 (Induced Subcomplex). There exists t± such that for all t > t+, Delp is an induced 
subcomplex of Del^t • 

Proof. For all t > 0, and each p G P, we have p G Vor mAp)- That is, input point are always their 
own nearest neighbors as weights increase. This implies that the vertices of Delp are all present in 
Del Mt • 

Suppose for contradiction that some simplex a is in Delp \ DeiM t for some t > r 2 , where r is 
the circumradius of a. Then there must be some Steiner point q encroaching the orthoball of a 
at time t. Let x be the circumcenter of a. Since a is composed only of input points, x is also the 
orthocenter of a for all times t. For all p G a, we have 

TT p (x, t) 2 = \\p — x\\ 2 — t = r 2 — t < < \\q — x\\ 2 = Tr q (x, t) 2 . 

It follows that q does not encroach on a at time t. This contradiction implies that Delp C DelM^ 
as long as is greater than the largest squared circumradius of Delp. 

To show that Delp is an induced subcomplex and complete the proof, we observe that because 
Delp and DelM t are embedded simplicial complexes covering the convex closure of P, there can be 
no simplices in Del^t containing only vertices of P that are not already included in Delp. □ 
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Finally, the last important fact to check is that there is not too much extra work to put the 
points in a bounding domain. In meshing, this slack between the bounding box and the input is 
called scaffolding |16j . 

Lemma 3 (Bounded Scaffolding). The number of simplices removed in the final step of the algo- 
rithm is 0(f). That is, \Bel Mi , \ Del P | = O(f). 

Proof. Any simplex a £ Deljvf* \Delp can be written as a disjoint union a = SUT where S £ ch(M) 
and T £ ch(P). By construction |ch(M)| is a constant. Since ch(P) C Delp, we have |ch(P)| = 
O(f). It follows that there can be at most |ch(M)| • |Delp| = O(f) such simplices. □ 

4.2 Running time analysis 

Theorem 4 (Running Time). Given n points P C M. d , the MeshVoronoi algorithm constructs 
the Voronoi diagram of P in 0(n log A log n) time. 

Before proceeding to the proof of the running time guarantee, we will first bound the num- 
ber of combinatorial changes of the weighted Voronoi diagram and thus also the number of heap 
operations. These are the main technical points of the analysis. 

We first prove a relatively standard packing bound on the number of Voronoi cells of Vor^ that 
can intersect the Voronoi cell of an input point. 

Lemma 5 (Packing). Let M be a t -well- spaced superset of P satisfying the feature size condition 
that fp(q) > K{m(q) for all q £ M for some constants r and K. Let p be any point of P. Then, 
the number of points q £ M such that VoiM(q) H Vorp(p) ^ ill is 0(log(aspectp(p))). 

Proof. The proof will be by a volume packing argument. Let M p denote the set of points q such 
that Vorjyf (<?) H Vorp(p) ^ 0. For any q £ M p and x £ Voi M(q) H Vorp(p), we derive the following 
bound on the distance between p and q. 

\\p — q\\ < \\q — x\\ + ||p — x\\ [by the triangle inequality] 

< \\q — x\\ + fp(x) \p is the nearest neighbor of x in P ] 

< 2\\q — x\\ + fp(q) [fp is 1-Lipscitz] 

< 2rfj\/((7) + fp(q) [since M is r- well-spaced] 

< (2t + K)¥m (q) [by the feature size condition] 

Define 7 to be | the constant in the above bound, i.e. 7 := 4j ^j ) _ 2T . For each integer i, define the 
spherical shell Si as follows. 

Si := {y £ R d I (1 + 7 r 1 < \\P ~ y\\ < (1 + 1)1- 

For each q £ M p , define the ball B q := ball(g, j\\p — q\\). The definition of B q and the bound on 
\\p-q\\ imply 

B q C baU(g, \hl{q)) C Vor M {q), 

and so the balls B q are pairwise disjoint. For q £ M p f)Ai, we further get that B q C ball(p, (l+7) l+1 ). 
Thus, 

Vol(ball(p,(l + 7 ) i+1 )) > Vol I |J Vol(B ff ) J > | ^|Vol(ball( 9 , 7(1 + 7) i_1 ))- 
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It follows that \Ai\ < (1 + 7) 2d /7 d . The max and min values of i such that Ai n M p is nonempty 
correspond to the nearest and farthest Voronoi neighbors of p in Vorp, so there can be only 
0(log(aspectp(p))) nonempty sets Ai n M p . This completes the proof as we have shown M p can be 
decomposed into 0(log(aspectp(p))) sets, each of constant size. □ 

Counting Flips The main challenge in the analysis is to bound the number of flips that happen 
in the transformation from the Voronoi diagram of M to the Voronoi diagram of P. The key to 
bounding this number is to observe that each such flip is witnessed by the intersection of a fc-face 
of Votm an d a (d — k)-face of Vorp for some k. Thus, we count these intersections instead. This 
intuition is made precise in the following lemmas. The first bounds the number of flips performed 
by the algorithm. The latter bounds the number of potential flips considered by the algorithm, 
since not all potential flips are performed. 

Lemma 6 (Flip Bound). Given n points P C M. d , the MeshVoronoi algorithm performs 0(f log A) 
flips, where f = |Vorp|. 

Proof. We observe that a flip occurs exactly when the weights cause some pair of adjacent simplices 
to encroach on each other. That is, the two simplices share a common orthoball. Let U be the 
d + 2 vertices comprising these simplices and let c be the center of their orthoball. The input points 
of U have weight 0, so they are equidistant from c. The weights of the Steiner points in U are 
all equal and thus these points are all also equidistant from c. So, the point c is the intersection 
of a fc-face of Vorp and a (d — /c)-face Vor^f, where M is the bounded aspect ratio superset of P 
constructed by the algorithm. So, it will suffice to bound the number of such intersections to get 
an upper bound on the number of flips. 

Suppose we have a /c-face of Vorp intersecting a (d — /c)-face Yotm- Then, it must be that there 
is a fc-face of Vorp intersecting a d-face of Vor^ . Since the d- faces of Vorjvf have only a constant 
number of faces, it will suffice to bound intersections between the d- faces of Vorp and the d- faces 
of Vorjvf- From Lemma [5j there are at most O(logA) such intersections. □ 

Lemma 7 (Potential Flip Bound). Given n points P C M. d , the MeshVoronoi algorithm sees 
0(/logA) potential flips, where f = |Vorp|. 

Proof. At the start of the algorithm, there is one potential flip for each facet of Vorjv/. This 
is at most 0{n log A). During the rest of the algorithm, there are at most ( rf ^ 2 ) = 0{d 2 ) new 
potential flips each time a real flip occurs, one for each new facet that appears. By Lemma [6j this 
is 0(f log A). □ 

The main result. We are now ready to prove running time guarantee, Theorem [4j 

Proof of Theorem^ It will suffice to bound the running time of each phase of the algorithm. The 
preprocessing takes 0(n log A) from the running time of Sparse Voronoi Refinement. Seeding the 
heap also takes 0(n log A) time as each facet is checked in constant time and the amortized cost of 
heap insertion is 0(1). There are at most two heap operations for each potential flip, one insertion 
and one deletemin. The total number of potential flips is 0(f log A) as shown in Lemma[7| Deleting 
the minimum element from a heap with 0(/logA) elements requires 0(log(n log A)) = O(logra) 
time. Thus, the total time for all heap operations is 0(n log Alogn) as desired. □ 
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5 Conclusion 



The algorithm we have presented is a direct combination of Delaunay mesh generation and kinetic 
data structures. The output-sensitive running time depends on the log of the spread. This is the 
usual cost of doing a geometric divide and conquer. It remains an interesting question if it is 
possible to replace the log A term with a logn by some more combinatorial divide and conquer. 
The other log-term coming from the heap operations may also permit some improvement as it is 
clear that many flips are geometrically independent and so their ordering is not strict. 

Another possible direction of future work is to exploit hierarchical meshes to keep the number 
of Steiner points linear [18]. However, it is not known how to leverage this into an improvement 
over the bounds presented here. At best it replaces the log A with max{logA,n}. 

Going forward, we hope to extend the methods here to the convex hull problem. Likely, this 
will require significant new ideas, but it may at least be possible to get significant improvements 
for convex hulls in the presence of bounded spread. 
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